Atomic-scale simulations of nanocrystalline metals 
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Nanocrystalline metals, i.e. metals in which the grain size is in the nanometer range, have a range 
of technologically interesting properties including increased hardness and yield strength. We present 
atomic-scale simulations of the plastic behavior of nanocrystalline copper. The simulations show that 
the main deformation mode is sliding in the grain boundaries through a large number of uncorrelated 
events, where a few atoms (or a few tens of atoms) slide with respect to each other. Little dislocation 
activity is seen in the grain interiors. The localization of the deformation to the grain boundaries leads 
to a hardening as the grain size is increased (reverse Hall-Petch effect), implying a maximum in hardness 
for a grain size above the ones studied here. We investigate the effects of varying temperature, strain rate 
and porosity, and discuss the relation to recent experiments. At increasing temperatures the material 
becomes softer in both the plastic and elastic regime. Porosity in the samples result in a softening of 
the material, this may be a significant effect in many experiments. 
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I. INTRODUCTION 

The modeling of the mechanical properties of every- 
day materials is a very challenging problem. The main 
difficulty is the vastly different length and time scales at 
which the various processes occur during deformation — 
ranging from the Angstrom and sub-picosecond scales of 
the atomic processes, to beyond the millimeter and sec- 
ond scales of the macroscopic deformation. Naturally, 
very different modeling techniques are required to model 
phenomena at so different scales. Atomic-scale simu- 
lations (typically molecular dynamics) can handle time 
scales of up to p|rfew nanoseconds and system sizes of 
up to 10® atomsoD, although one is typically limited to 
significantly smaller system sizes and simulation times 
by the available computer resources, and by the need to 
repeat the simulations at different conditions. 

At larger length and time scales, one is forced to aban- 
don the atomistic description of the material. One op- 
tion is dislocation dynamics, where the fundamental unit 
of the simulation is the dislocation.B'Q The fundamental 
idea of dislocation dynamics is similar to that of molec- 
ular dynamics: by calculating the forces on the dislo- 
cation (from the interactions with each other as well as 
from any applied stress), the equations of motions can be 
solved numerically. The modeling becomes significantly 
more complicated than molecular dynamics, mainly be- 
cause dislocations are lines whereas atoms are points. 
For that reason, dislocation dynamics is often done in 
two dimepsions, where the dislocations become point 
defectsJa'O One difficulty encountered when the atomic- 
scale description of the material is abandoned is how to 
treat inherently atomic-scale processes. In the case of 
dislocation dynamics, it will mainly be processes such as 



dislocation annihilation and core effects when the dislo- 
cations are intersecting. Parameters describing such pro- 
cesses must be extracted from experiments, or cajguiated 
directly using selected atomic-scale simulations. Hu^lI 

At still larger length and time scales, continuum mod- 
els must be used. This will typically be some kind of 
finite-element based calculation, where the plastic behav- 
ior is described by a constitutive relation. 

A number of advanced simulation techniques have been 
proposed, where simulation techniques on several length 
scales are combined into a single multiscale simulation. 
Various ways of combining finite element calculations of 
the long-range elastic fields with atonric-scale simulations 
of some ijeffioas of the material have appeared in the 
literatureJlHlla Unfortunately, most of these methods are 
limited to quasi-static sinrulations, i.e. to zero tempera- 
ture. 

Various simulation techniques have been proposed to 
extend the time scales beyond the nanosecond scale 
reachable by molecular dynamics. Sonre nrethods are 
based on "accelerating" the nricroscopic dynamics^of the 
system by modifying the potential energy surfacdi^ or by 
extrapolating from higher to lower temperatures .ILl 
methods, such as the "nudged elastic tf J|i4" methocK 
and constrained molecular dynamics,Ej £3 allow one to 
determine energy barriers, from which reaction rates of 
slow processes can be determined. 

A different way of addressing the length-scale challenge 
is by simulating suitably chosen problems, where the 
characteristic length scale is within the reach of atomic- 
scale modeling. This has the advantage that the entire 
problem can be addressed at the atomic scale, remov- 
ing the need for the a priori assumptions about which 
processes should be included in the simulation — some- 



times unexpected atomic-scale processes are seen in such 
simulations, processes that could not have been included 
in a coarser-scaled simulation as its importance was not 
suspected. 

One such atomic-scale problem is the mechanical defor- 
mation of nanocrystalline metals, i.e. metals where the 
grain size is in the nanometer range. Nanocrystalline 
metals have recently received much interest because they 
may have mechanical, chemical and physical properties 
different from their coarse-grained counterparts. For ex- 
ample, the hardness and yield stress may increase 5-10 
times when the grain sizeis reduced from the macroscopic 
to the nanometer rangeBjcJ 

Recently, computer simulations of the structure 
of nanocry^alline metals— asd semiconductors, and of 
their elasticEZl and plasticEilS properties have appeared 
in the literature. In previous papers, we described the 
plastic deformation of nanocrystalline copper at zero 
temperatureoo In this paper we focus on the elastic 
and plastic properties of nanocrystalline metals, in par- 
ticular copper, at finite temperature. We find that the 
materials have a very high yield stress, and that the yield 
stress decrease with decreasing grain size (reverse Hall- 
Petch effect). The main deformation mode is found to 
be localized sliding in the grain boundaries. 

The high yield stress and hardness of nanocrys- 
tallinft-metals is generally attributed to the Hall-Petch 
effect ,HE3 where the hardness increases with the inverse 
square root of the grain size. The Hall-Petch effect is 
generally assumed to be caused by the grain boundaries 
acting as barriers to the dislocation motion, thus hard- 
ening the material. The detailed mechanism behind this 
behavior is still under debate.Ell A cessation or reversal of 
the Hall-Petch effect will therefore limit the hardness and 
strength that can be obtained in nanocrystalline metals 
by further refining of the grain size. There are a num- 
ber of observations of a reverse Hall-Petchj-afEect, i.e. of 
a softening when the grain size is reduced. EaliU The in- 
terpretation of these results have generated some con- 
troversy. It is at present not clear if the experimentally 
reported reverse Hall-Petch effect is an intrinsic effect or 
if it is caused by reduced sample quality at the finest 
grain sizes. The computer simulations presented here 
show that an intrinsic effect is clearly possible. 

The structure of the paper is as follows. In section || we 
discuss the setup of the nanocrystalline model systems. 
Section HI discusses the simulation and analysis methods 
used. The simulation results are presented in section PV| , 
and subsequently discussed in section 0. 



II. SIMULATION SETUP 

In order to obtain realistic results in our simulations, 
and to be able to compare our simulations with the avail- 
able experimental data, we have attempted to produce 
systems with realistic grain structures. Unfortunately, 



the microscopic structure is not very well characterized 
experimentally, and depends on the way the nanocrys- 
talline metal was prepared. We have tried to create sys- 
tems that mimic what is known about the grain structure 
of nanocrystalline metals generated by inert gas conden- 
sation. The grains seem to be essentially equiaxed, sep- 
arated by narrow grain boundaries. The grains are es- 
sentially dislojcaiipn free. The grain size distribution is 
log-normal .t^OE^I 



A. Construction of the initial configuration 

In our simulationS|-t|he grains are produced using a 
Voronoi construction:Ej a set of grain centers are chosen 
at random, and the part of space closer to a given center 
than to any other center is filled with atoms in a ran- 
domly rotated face-centered cubic (fee) lattice. Periodic 
boundary conditions are imposed on the computational 
cell. This procedure generates systems without texture 
and with random grain boundaries. Effects of texture 
could easily be included by introducing preferred orien- 
tations of the grains. In the limit of a large number of 
grains, the Voronoi construction will generate a-gpain size 
distribution close to a log-normal distributional 

In the grain boundaries thus generated, it is possible 
that two atoms from two different grains get too close to 
each other. In such cases one of the atoms is removed to 
prevent unphysically large energies and forces as the sim- 
ulation is started. To obtain more relaxed grain bound- 
aries the system is annealed for 10 000 timesteps (50 ps) 
at 300 K, followed by an energy minimization. This pro- 
cedure is important to allow unfavorable local atomic 
configurations to relax. 

To investigate whether the parameters of the annealing 
procedure are critical, we have annealed the same sys- 
tem for 50 and 100 ps at 300 K, and for 50 ps at 600 K. 
We have compared the mechanical properties of these 
systems with those of an identical system without an- 
nealing. We find that the annealing is important (the 
unannealed system was softer), but the parameters of 
the annealing are not important within the parameter 
space investigated. 

A aimilar generation procedure has been repoijted by 
ChenjHj by D'Agostino and Van Sswgenhoven ,c3 and 
by Van Swygenhoven and Caro.LJLS A different|-ap-| 
proach was proposed by Phillpot, Wolf and Gleiter£3'E2l 
a nanocrystalline metal is generated by a computer sim- 
ulation where a liquid is solidified in the presence of crys- 
tal nuclei, i.e. small spheres of atoms held fixed in crys- 
talline positions. The system was then quenched, and 
the liquid crystallized around the seeds, thus creating a 
nanocrystalline metal. In the reported simulations, the 
positions and orientations of the seeds were determinis- 
tically chosen to produce eight grains of equal size and 
with known grain boundaries, but the method can natu- 
rally be modified to allow randomly placed and oriented 



seeds. The main drawback of this procedure is the large 
number of defects (mainly stacking faults) introduced in 
the grains by the rapid solidification. The stacking faults 
are clearly seen in the resulting nanocrystalline metal 
(Fig. 7 of Ref . E^) . The appearance of a large number of 
stacking faults was also seen in the solidification of large 
clusters even if the cooling-is done as slowly as possible 
in atomistic simulations.OO 



B. Structures 

A typical system (after annealing and cooling to zero 
Kelvin) with a grain size of 5.2 nm is shown in Fig. n. 
The atoms have been color coded according to the local 
crystal structure, as determined by the Common Neigh- 
bor Analysis (see section niC| ). In Figure g the radial 
distribution function (RDBF) g(r) for the same system is 
shown. It is defined as the average number of atoms per 
volume at the distance r from a given atom. The RDBF 
is seen to differ from that of a perfect fee crystal in two 
ways. First, the peaks are not sharp delta functions, but 
are broadened somewhat. This broadening is in part due 
to strain fields in the grains (probably originating from 
the grain boundaries), and in part due to atoms in or 
near the grain boundaries sitting close to (but not at) 
the lattice positions. The second difference is seen in the 
inset: the RDBF does not go to zero between the peaks. 
It is the signature of some disorder which in this case 
comes from the grain boundaries. 

Experimentally, information about the RDBF can be 
obtained from X-ray absorption fine structure (XAFS) 
measurements. This method has been used to mea- 
sure the average coordination number of Cu atoms 
in nanocrystalline Cu, finding coordination numbers of 
11.8±0.3 and 11.9±0.3JJi samples with 34 nm and 13 nm 
grain size, respectively.c3 From these results the average 
coordination number of the atoms in the grain boundaries 
was estimated to 11.4 ± 1.2, i.e. within the experimental 
uncertainty it is the same as in the bulk. Integrating the 
first peak of the calculated RDBF (Fig. ||) gives an av- 
erage coordination number of 11.9 ± 0.15. As the RDBF 
does not go to zero between the first two peaks, it is 
not clear where the upper limit of the integration should 
be chosen, hence the uncertainty. The value given is for 
an upper limit of 3.125 A. There is thus excellent (but 
perhaps rather trivial) agreement between the calculated 
and the experimental coordination numbers. 

Numerical studies have shown that the Voronoi con- 
struction results in a grain size distribution that it is 
well described by a log-normal distribution (although for 
more than 5000 grains aj-two-parameter gamma distri- 
bution gives a better fit).L3 In Fig. we show the grain 
size distributions in our simulations with intended av- 
erage grain sizes of 3.28 and 5.21 nanometers. The ob- 
served distributions are consistent with log-normal distri- 
butions, although due to the rather low number of grains 



it is not possible to distinguish between a log-normal or 
a normal distribution. 



III. SIMULATION METHODS 

We model the interactions between the atoms us- 
mg a maap|body Effective Medium Theory (EMT) 



potential. e3'L1I EMT gives a realistic description of the 
metallic bonding, in particular in fee metals and their 
alloys. Computationally, it is not much more demanding 
than pair potentials, but due to its many-body nature 
it gives a far more realistic description of the materials 
properties. 

The systems can be deformed by rescaling the coordi- 
nates along a direction in space (in the following referred 
to as the z direction). During this deformation either 
a conventional molecular dynamics (MD) algorithm or 
a minimization algorithm is used to update the atomic 
positions in response to the deformation. 

In molecular dynamics, the timestep used when inte- 
grating the equation of motion must be short compared 
to the typical phonon frequencies in the system. We use 
a timestep of 5 fs, safely below the value where the dy- 
namics becomes unstable. A consequence of the short 
timesteps involved in MD simulations is that only brief 
periods of time can be simulated. For the size of systems 
discussed here 1 ns (200000 timesteps) is for all practical 
purposes an upper limit, although for repeated simula- 
tions 0.1 ns is a more realistic limit. 

One consequence of the short time-scale is that very 
high strain rates are required to get any reasonable de- 
formation within the available time, a typical strain rate 
in the simulations reported here is 5 x lO^s^^. This is 
very high, but still the ends of the system separates at 
velocities far below the speed of sound. We have inves- 
tigated the effects of varying the strain rate, see section 
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As another consequence of the short time-scale, slower 
processes will not be seen in the simulations. In par- 
ticular, most diffusional processes will be unobservable 
(this would also be the case in experiments performed 
at these high strain rates). However, measurements of 
diffusional creep (Coble creep) in nanocrystalline iH«te|ls 
indicate that diffusional creep is not a large effect E3H 

Some of the simulations were performed using a mini- 
mization procedure, i.e. the system is kept in a local en- 
ergy minimum while it is deformed. In such a simulation 
time is not defined, since we are not solving an equation 
of motion. In a similar way, time should not be relevant 
in an experiment performed truly at zero temperature, 
since there will be no thermally activated processes, and 
thus no way the sample can leave a local energy mini- 
mum (neglecting quantum tunelling). So the strain rate 
will not have an effect on the results, providing it is low 
enough to prevent a heating of the sample, and provid- 
ing the minimization procedure is fully converged. The 



minimization simulations can thus be seen as a model 
for idealized experiments at zero temperature in the low 
strain rate limit, where there is time for the heat gener- 
ated by the deformation to be removed. 



A. The minimization procedure 

To simulate deformation at zero temperature a min- 
imization procedure was used to keep the system in or 
near a local minimum in energy at all times. The defor- 
mation and minimization was done simultaneously. The 
minimization algorithm is a modified molecular dynam- 
ics simulation.Ej After each MD timcstcp the dot product 
between the momentum and the force is calculated for 
each atom. Any atom where the dot product is negative 
gets its momentum zeroed, as it is moving in a direction 
where the potential energy is increasing. In this way the 
kinetic energy of an atom is removed when the potential 
energy is close to a local minimum along the direction 
the atom is moving. This minimization procedure quickly 
brings the system close to a local minimum in energy, but 
a full convergence is not obtained, as it would require a 
number of timesteps at least as great as the number of 
degrees of freedom in the system. However, we find only 
little change in the development of the system, when we 
increase the number of minimization steps. 

At each timestep the system is deformed by a tiny scal- 
ing of the coordinates, the z coordinates are multiplied 
by 1 -t- e, the x and y coordinates by 1 — i^e, where e 
is a very small ruimber, chosen to produce the desired 
deformation rateJ£j The constant v is an "approximate 
Poisson's ratio". A Monte Carlo algorithm is used to 
optimize the two lateral dimensions of the system: after 
every 20th timestep a change in the lateral dimensions 
is proposed. If the change result in a reduction of the 
energy it is accepted, otherwise it is discarded. In this 
way the exact value chosen for v becomes uncritical, as 
the Monte Carlo algorithm governs the contraction in the 
lateral directions. The use of ;/ ^ is just for computa- 
tional efhciency. We used v = 0.4 as this was between the 
optimal value in the part of the simulation where elas- 
tic deformation dominates (Poisson's ratio « 0.3) and in 
the part where plastic deformation dominates (y « 0.5 
as volume is then conserved). 

A few simulations were repeated using a Conjugate 
Gradient (CG) minimization instead of the MD mini- 
mization algorithm. The two algorithms were approx- 
imately equally efficient in these simulations, provided 
that the CG algorithm was restarted approximately ev- 
ery 20 line minimizations. Otherwise the CG algorithm 
will not minimize twice along the same direction in the 
SA^-dimensional configuration space. 



B. Molecular dynamics at finite temperature 

At finite temperatures a conventional molecular dy- 
namics algorithm is used, where the Newtonian equations 
of motion for the atoms are solved numerically. During 
the simulation the computational box was stretched as 
described above. The Monte Carlo algorithm optimiz- 
ing the lateral dimensions is the conventional Metropolis 
algorithm. 

Before the deformation is applied the system is heated 
to the desired temperature by a short Moleaular Dynam- 
ics simulations using Langevin dynamics ,E3 i.e. where a 
friction and a fluctuating force is added to the equation of 
motion of the atoms. When the desired temperature has 
been reached (after approximately 10 ps), the simyilations 
is performed using the velocity Verlet algorithm.EZI Dur- 
ing the deformation process the internal energy increases 
by the work performed on the system. This amounts in 
practice to a small heating of the system of the order of 
-30 K. 



C. Analysis of the results 

While the simulation is performed, the stress field is 
regularly computed. The stress tensor is the derivative 
of the free energy of the system with respect to the strain. 
The effective medium theory allows us to define an energy 
per atom, which allows us to define an "atomic" stress 
for each atom. The stress is a suitable derivative-pt-the 
energy with respect to the interatomic distancesr* 
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where ai^ap is the a, (3 component of the stress tensor for 
atom z, Vi is the volume assigned to atom i {^^ Vi = V, 
where V is the total volume of the system), rrii is the 
mass of atom ?, pi^a is the a component of its momentum 
and Tij is the distance between atoms i and j (r^-.Q is a 
component of the vector from atom i to j). 

The atomic stress tensor cannot be uniquely defined. 
Eq. is based on JJne_yirial Theorem, but other defini- 
tions are possible .|£J"E3 When the atomic stress is aver- 
aged over a region of space the various jdafijiitions quickly 
converge to a macroscopic stress field.ESS 

During the simulation, stress-strain curves are calcu- 
lated by averaging the atomic stresses over the entire 
system. 

To facilitate the analysis of the simulations the local 
atomic order was examined using an .aljarithm known 
as Common Neighbor Analysis (CNA).E3^lj In this algo- 
rithm the bonds between an atom and its neighbors are 
examined to determine the crystal structure (two atoms 
are considered to be bonded if they are closer together 
than a cutoff distance chosen between the first two peaks 
in the radial distribution function). Bonds are classified 



by three integers (ijk). The first integer i is the number 
of common neighbors, i.e. atoms bonded to both atoms 
in the bond under consideration. The second integer j 
is the number of bonds between these common neigh- 
bors. The third integer k is the longest chain that can 
be formed by these bonds. 

The number and type (ijk) of bonds that an atom 
has determines the local crystal structure. For example, 
atoms in a perfect fee crystal have 12 bonds of type 421, 
whereas atoms in a perfect hep crystal have six bonds 
of type 421 and six of type 422. In these simulations 
we have mainly used CNA to classify atoms into three 
classes: fee, hep and "other" , where the "other" class is 
atoms that have a different number of bonds than 12, or 
that have bonds that are not of type 421 or 422. 

The use of CNA makes dislocations, grain boundaries 
and stacking faults visible in the simulation. Intrinsic 
stacking faults appear as two adjacent 111 planes of hep 
atoms, extrinsic stacking faults are two 111 planes of 
hep atoms separated by a single 111 plane of fee atoms, 
whereas twin boundaries will be seen as a single 111 
plane of hep atoms. Dislocation cores and grain bound- 
aries consist of atoms in the "other" class, although grain 
boundaries also contain a low number of hep atoms. 

When analyzing simulations made at finite tempera- 
tures, the lattice vibrations may interfere with the CNA. 
For simulations at temperatures up to 300 K it is suffi- 
cient to choose the cutoff distance carefully, but for higher 
temperatures it is necessary to precede the CNA analysis 
by a short minimization (short enough to remove most 
of the lattice vibrations, but not long enough to allow 
dislocations etc. to move). For consistency, such a min- 
imization procedure was used for all finite-temperature 
simulations regardless of temperature. 



IV. RESULTS 

In this section we first discuss the results of the sim- 
ulations based on the minimization algorithm (i.e. the 
zero-temperature results), and then we discuss the re- 
sults obtained using molecular dynamics at finite tem- 
peratures. 



A. Simulations at zero temperature 

During the deformation, we calculate the average stress 
in the system as a function of the strain. Fig. ^ shows the 
obtained stress-strain curves from simulations at OK. We 
see a linear elastic region followed by a plastic region with 
almost constant stress. Similar results are found for pal- 
ladium (Fig. W. Each stress-strain curve shown in Figs. H 
and ^ are obtained by averaging over a number of simula- 
tions with different (randomly produced) grain structures 
with the same average grain size. A set of stress-strain 
curves from individual simulations are shown in Fig. 0. 



One of the rationales using a minimization procedure 
to study the deformations was the hope that the system 
would evolve through a series of local energy minima, sep- 
arated by discrete events when the applied deformation 
causes the minima to disappear. In this way, the simula- 
tion would have resulted in a unique deformation history 
for any given sample. However, the deformation turned 
out to happen through a very large number of very small 
processes, that could not be individually resolved by this 
procedure (see below). One symptom of this is, that 
the individual curves in Fig. @ are not completely repro- 
ducible. Any even minor change in the minimization pro- 
cedure, or a perturbation of the atomic coordinates, will 
result in a slightly different path through configuration 
space, and in different fluctuations in the stress-strain 
curves. Those differences are suppressed when average 
stress-strain curves are calculated, as in Fig. ]^ and would 
also disappear as the system size (and thus the number 
of grains) increase. 



1. Young's Modulus 

In the linear elastic region the Young's modulus is 
found to be around 90-105 GPa and it is increasing with 
increasing grain size. The experimentalrAtalue for macro- 
crystalline copper is 124 GPa at 300 kJ^ and the value 
found for single crystals using this potential is 150 GPa 
at K (Hill average calculated from the anisotropic elas- 
tic constants Cn = 173 GPa, C12 = 116 GPa, C44 — 
91 GPa). A similar reduction of Young's modulus is seen 
in simulations—pf nanocrystalline metals grown from a 
molten phase. tZi The low value is due to the large volurpsj 
fraction of the atoms being in the grain boundaries .eZI 
These atoms experience a different atomic environment, 
which could result in a reduction of the elastic moduli 
similar to what is seen in amorphous metals. This local 
reduction of the elastic constants in_grain boundaries is 
confirmed by atomistic simulations .L2 

Experimental measurements of the Young's modulus 
of high-quality (i.e. low-porosity) samples of nanocrys- 
talline copper and palladium show a reduction in Young's 
modulus of at mostj-a. few percent when correcting for the 
remaining porosity.Ej These results were obtained for sig- 
nificantly larger grain sizes than were used in the simula- 
tions. The reduction of Young's modulus that we observe 
in these simulation, will be difficult to detect experimen- 
tally due to the much lower volume fraction of atoms in 
the grain boundaries for typical grain sizes in high-quality 
samples (> 20 nm). 



2. Yield and flow stress 

The onset of plastic deformation is usually described 
by the yield stress ay, traditionally defined as the stress 



where the strain is 0.002 larger than what would be ex- 
pected from extrapolation from the elastic region. In 
these simulations the stress continues to increase after the 
yield point has been reached, until it reaches a plateau 
and becomes constant (or slightly decreasing). We call 
the level of the plateau the flow stress. 

Fig. 1^ shows the dependence of the yield and flow stress 
on the grain size. A clear reverse Hall-Petch effect is 
observed, i.e. a softening of the material as tJae grain size 
is reduced, as discussed in a previous papero 



3. Structural changes 

Fig. g shows the same system as Fig. |l|, but after 10 % 
deformation. Some stacking faults have appeared in the 
grains, they are caused by partial dislocations (Shockley 
partials) nucleating at the grain boundaries and moving 
through the grains. One such dislocation is seen in the 
figure. 

The radial distribution function (Fig. ^ has been 
changed somewhat by the deformation. The peaks have 
been broadened, this is mainly caused by the anisotropic 
stress fields in the sample. The "background level" be- 
tween the first two peaks has increased a little, indicat- 
ing a larger amount of disorder in the system. Increased 
disorder is also seen in Fig. H, where the grain bound- 
aries appear to have become slightly thicker compared 
to the initial configuration. This is confirmed by Fig. nG, 
showing the number of atoms in different local configura- 
tions before and after the deformation. We see how more 
atoms are neither fee nor hep after the deformation than 
before. 

A strong increase in the number of atoms near stacking 
faults (atoms in hep symmetry) is also seen in Figs. 
and 10, The stacking faults appear as partial dislocations 
move through the system, and they are thus the signature 
of dislocation activity. At zero temperature, we do not 
observe cases where a second partial dislocation erases 
the stacking faults (we observe only a very few atoms 
changing from local hep order to local fee order). We 
can therefore use the total number of hep-ordered atoms 
to estimate an upper bound on the amount of plastic 
deformation caused by the dislocations. 

If a dislocation with Burgers vector b runs through the 
entire system, the dimensions of the system are changed 
by b and the strain Ezz is thus bzL~^, where bz is the z 
component of the Burgers vector and Lz is the dimen- 
sion of the system in the z direction. If the dislocation 
only passes through a part of the system, the result- 
ing deformation is reduced by Acos(j){LxLy)~^ , where 
A is the area of the slip, L^ and Ly are the lateral di- 
mensions of the system and (j) is the angle between the 
slip plane and the xy plane. The contribution from a 
slip plane to the zz component of the strain is thus 
Szz = {bzL-^)Acos(j){LxLyy^ . 

The maximal value of bz is bsincf), where b — \b\, since 



b lies in the slip plane. The maximal strain from the 
slip thus becomes e^ax = bA{2V)~^ = aA{2V6V)~^ for 
(j) = 7r/4, as the Burgers vector of a Shockley partial 
is 6 = a/V6, where a is the lattice constant and V is 
the volume of the simulation cell. A slip plane of area A 
results in two {111} planes of hep atoms, i.e. 4A(\/3a^)^^ 
atoms. The total system contains 4:Va~^ atoms, so the 
fraction of hep atoms becomes n = Aa{\/6V)^^ . Hence 
a fraction n of hep atoms can at most have resulted in 
a strain of Emax = 2~^/^n. As the simulation generate 
at most 9% hep atoms during 10% deformation, we get 
that Emax ~ 3%, provided that all slip planes and Burgers 
vectors are ideally aligned. It is therefore clear that the 
main deformation mode is not by dislocation motion. 

Fig. ITl| illustrates how the main part of deformation 
has taken place. The atoms are colored according to 
their motion relative to the global stretching of the sys- 
tem. We clearly see that the upper parts of the grains 
have moved down and the lower parts up, relative to what 
would be expected in a homogeneous deformation. This 
shows that the grains do not stretch as much as in a ho- 
mogeneous deformation. On the other hand, it is seen 
that significant deformation has happened in the grain 
boundaries, as the atoms typically are moving up on one 
side and down on the other side of a grain boundary. An 
analysis of that deformation shows that it is in the form 
of a large number of apparently uncorrelated small slip- 
ping events, where a few atoms-Lor a few tens of atoms) 
move relatively to reach other,E3 i.e. not in the form of 
collective motion in the grain boundaries. A minor part 
of the plastic deformation is in the form of dislocation 
motion inside the grains. The slip planes are clearly seen 
in Fig. O, in particular in the large grain in the upper 
left part of the figure, where two dislocations have moved 
through the grain, and a third is on its way near one of 
the previous slip planes. 



B. Simulations at finite temperature 

The simulations were repeated with the same systems 
(i.e. the same initial grain structure) at finite temper- 
atures. Fig. |l2| shows the stress-strain curves for the 
same system at different temperatures. We clearly see 
a softening with increasing temperature. Rather large 
fluctuations are seen in the curves. These are mainly 
thermal fluctuations and fluctuations due to single major 
"events" in the systems (e.g. the nucleation of a dislo- 
cation). These fluctuations are only visible due to the 
small system size. 

The softening of the material with decreasing grain size 
is also observed in simulations at 300 K, see Fig. nj. 



1. Young's Modulus 

Young's Modulus {E) is the slope of the stress-strain 
curve in the linear region. When calculating Young's 
modulus from the simulation data, a compromise must be 
made between getting enough data point for a reliable fit, 
and staying within the clearly linear region. For the zero- 
temperature simulations, fitting Young's modulus to the 
data point for e < 0.1% satisfies both conditions, but for 
the finite-temperature simulations more data points are 
required. We have chosen to use data in the interval e < 
1%, this ensures that we have enough data for a reliable 
fit, but results in a slight underestimate of the Young's 
modulus, as some plastic deformation is beginning in this 
interval. The results of this procedure is shown in Fig. |l4| , 
showing Young's modulus for a single system with grain 
size d = 5.2 nm, simulated at different temperatures. For 
consistency, the larger strain interval has been used even 
for the T — simulation. Using the smaller interval 
(e < 0.1%) would result in E ~ 119 Gpa instead of E = 
100 GPa. 

The observed temperature dependence of E is approx- 
imately -72 MPa/K, which is somewhat larger than what 
has been observed experimeotally (-40 MPa/K) in copper 
with a grain size of 200nmJlJ This may be because the 
Young's modulus of the grain boundaries is more tem- 
perature sensitive than in the bulk, or it may be due to 
increased creep in the higher-temperature simulations, 
see the discussion in section [V B 4 . 



2. Yield and flow stress 

The Yield stress is again determined as the stress where 
the strain is 0.2% above what would be expected from ex- 
trapolation from the linear regime. The difhculties lead- 
ing to an underestimate of the Young's Modulus thus 
leads to an overestimate of the yield stress. The values 
for the yield stress obtained for the simulations at 300 K 
can therefore not be compared directly with the values 
obtained at OK, but values obtained at different grain 
sizes can of course be compared, as the same method 
was used to estimate the yield stress in all cases. The 
flow stress is a far more well-defined quantity, and direct 
comparison is possible between simulations at different 
temperatures. The variation of the yield and flow stresses 
with temperature is seen in Fig. HS. 

Fig. 16 shows the variation of both the yield- and the 
flow stress with grain size. As in the simulations at K 
a reverse Hall-Petch relationship is found. 



3. Structural changes 

The main deformation mode appears to be the same 
at zero and at finite temperatures. Figure |^ shows the 



atomic displacements, again the majority of the defor- 
mation has taken place in the grain boundaries, only a 
few slip planes are seen. 

The grain boundaries do not appear to increase as 
much in thickness as they to at OK. Fig. n^ shows the 
change in the fraction of atoms in different local environ- 
ments during the deformation of the same system at K 
and at 300 K. In both cases we clearly see an increase in 
the number of hep atoms (stacking faults) due to the mo- 
tion of dislocations through the grains, but the number 
of atoms in the grain boundaries increases significantly 
more at K than at 300 K (the increase is approximately 
twice as big at OK as at 300 K). The increase appears to 
be caused by the deformation in the grain boundaries. 
Apparently, the local disorder introduced in this way is 
partially annealed out at 300 K. 



4- Strain rate 

The finite-temperature simulations presented in this 
paper were performed at a strain rate of e = 5 x 10*s~^, 
unless otherwise is mentioned. 

In order to investigate the influence of the strain rate, 
we simulated the same deformation of the same system 
using different strain rates in the range 2.5 x 10^s~^ - 
1.0 X 10^°s^^, see Fig. |l9|. A strong dependence on the 
strain rate is seen for strain rates above 1 x 10^s~^ . Below 
this "critical strain rate" the strain rate dependence on 
the stress-strain curves is far less pronounced. Fig. EG 
confirms this impression. It shows the yield and flow 
stress as a function of the strain rate. Experiments on 
ultrafine-grained (d w 300 nm) Cu and Ni show a clear 
strain rate dependence on the yield and flow stresses at 
high strain ratesO 

Perhaps surprisingly, the Young's modulus appears to 
depend on the strain rate as well (Fig. 11^). This indi- 
cates that some kind of plastic deformation occurs in the 
"linear elastic" region. This is confirmed by stopping a 
simulation while the system still appears to be in the 
elastic region, and then allowing it to contract until the 
stresses are zero. The system does not regain the original 
length: plastic deformation has occurred. 

To examine the time-scale over which this deformation 
occurs, a configuration was extracted from the simulation 
at e = 2.5 x 10^s~^ after 0.4% deformation. The system 
was held at a fixed length for 300 ps while the stress is 
monitored, see Fig. ^. The stress is seen to decrease with 
a characteristic time of approximately 100 ps. By plot- 
ting the atomic motion in a plot similar to Figs. ^ and 
pT) , it is seen that the relaxation is due to small amounts 
of plastic deformation in the grain boundaries. The con- 
sequence of this is that the systems do not have time to 
relax completely during the simulations, explaining the 
observed strain rate dependence. In order to allow for 
complete relaxation of the systems, strain rates far be- 
low what is practically possible with MD simulations are 



required. 



5. Grain rotation 



Grain rotation has previously|-heen reported in simu- 
lations of nanocrystalline nickel.tS We have investigated 
the rotation of the grains during some of the simulations, 
the results are summarized in Fig. E3. The figure shows 
the rotation of five randomly selected grains as a function 
of strain and temperature. The rotations were identified 
by a three-dimensional Fourier transform of the positions 
of the atoms in the grains. 

We see that the grain rotation increases with increas- 
ing temperature. There is a large variation between how 
much the individual grains rotate. The grains with the 
largest rotations keep the same axis of rotation during 
the entire deformation, whereas the grains that only ro- 
tate a little have a varying axis of rotation. Probably 
some grains are in a local environment where a signifi- 
cant rotation results in an advantageous deformation of 
the sample which reduces the stress. Other grains are 
randomly rotated as the many small deformation pro- 
cesses in the grain boundaries occur. 



6. Porosity 

As the observed reverse Hall-Fetch effect is often ex- 
plained as an artifact of sample porosity (see section 
VA), we found it relevant to study how pores influence 



the mechanical properties. The void structure in exper- 
imentally produced samples is usually not well known, 
so we chose to study several different types of voids. In 
all cases the voids resulted in a reduction of both the 
Young's modulus and the flow stress, see Fig. g^. 

Elliptical voids. These crack-like voids were created 
by removing all atoms within an oblate ellipsoid with 
an axis ratio of 3.16. The short axis can be oriented 
along the pulling direction (the z-axis) or perpendicular 
to it (the X-axis) . The former orientation corresponds to 
cracks that are activated by the applied stress field, the 
effect of these cracks is therefore expected to be much 
larger than the effect of the "inactive" cracks. This is 



clearly seen in Fig. |23 
Missing grains 



There have been reports of pore sizgs 
comparable to (and proportional to) the grain size.E3'LHI 
To emulate this, we have tried to remove whole grains 
from the system. As the grains are approximately 
equiaxed, it is not surprising that the effect of removing 
a grain is intermediate between the effects of removing 
ellipsoids in the two orientations, provided that approxi- 
mately the same number of atoms are removed. 

Missing grain boundary atoms. In samples experi- 
mentally produced by compacting a powder, it is rea- 
sonable to assume that the porosity will mainly be in 
the form of (possibly gas-filled) voids between the grains. 



There is also *3jp§ experimental evidence that this is in- 
deed the case.L30 To emulate this, we have removed all 
atoms in the grain boundaries within one or nine spheri- 
cal regions in the sample, creating one large or nine small 
voids in the grain boundaries. This type of voids have 
the largest effect on the materials properties, giving a 
reduction of 35-40% in the Young's modulus and flow 
stress for a 12.5% porosity. It seems rather natural that 
a large effect is obtained with the voids concentrated in 
the grain boundaries since we know that the main part 
of the deformation is carried by these boundaries. 



V. DISCUSSION 

A. The reverse Hall-Petch effect 

A reverse Hall-Petch effect in nanocrystalline copper 
was first observeiip nanocrystalline Cu and Pd by Chok- 
shi et al. in 1989.E2I Since then, there have been nuroerpus 
observations of softening at very small grain sizes. EJll^ 

The reverse Hall-Petch effect seems to depend strongly 
on the sample preparation technique used and on the 
sample history, perhaps indicating that in most cases the 
reverse Hall-Petch effect is caused by various kinds of 
defects in the samples. Surface defects alone have been 
shown to be able to decrease thepSj|Uength of nanocrys- 
talline metals by a factor of five|£jO and recent stud- 
ies have shown that even very small amounts of. poros- 
ity can have a dramatic effect on the stpepgthEj'Ej Im- 
proved inert gas condensation techniqueal3 have reduced 
the porosity resulting in samples with densities above 
98% of the fully-dense value. In these samples the ordi- 
nary Hall-Petch effeat is seen to continue down to grain 
sizes around 15 nm.E3 There are only few data points be- 
low that grain size, but apparently no further increase 
in the hardness is seen. It is suggested that most of the 
observations of a reverse Hall-Petch effect in naxiocrys- 
tallinc metals are a result of poor sample qwaJity.LJ This 
impression is supported by literature studieg230 indicat- 
ing that the reverse Hall-Petch effect is mainly seen when 
the grain size is varied by repeated annealing of a single 
sample, whereas an ordinary Hall-Petch relationship is 
seen when as-prepared samples are used. 

However, there does seem to be a deviation from the 
Hall-Petch effect for grain sizes below approximately 
15 nm, where the Hall-Petch slope is seen to decrease 
or vanish in samples produced with various techniques. 
This is seen in Cu samples produced by inert gas con- 
densation followed by warm compaction (sample densi- 
ties above 98%)E3 and in electroplated Ni (claimed to be 
porosity free).L3 

There are theoretical arguments for expecting that the 
Hall-Petch relation ceases to be valid for grain sizes below 
^20 nm: as the grain size becomes too small, dislocation 
pileups are no longer possible, and the usualpexplanation 
for Hall-Petch behavior does not applySoEEl 



Many models have been proposed to explain why a 
revsrse Hall-Petch effect is sometimes seen. Chokshi et 
alx3 proposed that enhanced Coble creep, i.e. creep by 
diffusion in the grain boundaries, should result in a soft- 
ening at the smallest grain sizes as the creep rate in- 
creases with decreasing grain size (d) as d^^. Direct 
measuxements of the creep rate have however ruled this 

out.E3e 

It has been suggested that the grain boundaries in 
nanocrystalline metals have a different structure, making 
them more transmrpBt to dislocations than "ordinary" 
grain boundaries .ESla If it becomes possible for the dis- 
locations to run through several grains as the grain size is 
reduced, the Hall-Petch relations would break down. In 
our simulations, we have not observed dislocations mov- 
ing through more than one grain. 

If the Hall-Petch effect is explained by appealing to dis- 
location sources in the grain boundaries, the Hall-Petch 
relationship is expected to break down when the grain 
sizes becomes so low that there jUxe no longer dislocation 
sources in all grain boundariesEj (assuming a constant 
density of disioiftation sources in the grain boundaries). 

Hahn et aLES'lZl suggest that the reverse Hall-Petch ef- 
fect is caused by deformation in the grain boundaries. If 
a grain boundary slides, stress concentrations build up 
where the grain boundary ends, limiting further sliding. 
Substantial sliding on a macroscopic scale occurs when 
sliding occurs on slide planes consisting of many aligned 
grain boundaries; the sliding is hindered by the rough- 
ness of the slide plane due to its consisting of many grain 
boundaries. As the grain size is reduced and becomes 
comparable to the grain boundary width, the roughness 
of such slide planes decreases and the stress required 
for mesoscopic sliding decreases. This would result in 
a reverse Hall-Petch effect. They estimate the transition 
from normal to reverse Hall-Petch effect to occur at grain 
sizes near 50 nm for Cu.e3 

The simulations reported in the present paper indicate 
that the main deformation mechanism at these grain sizes 
is indeed sliding in the grain boundaries. However, it is 
not clear if the proposed "collective" sliding events are 
occurring, it appears that sliding occurs on individual 
grain boundaries, and that the resulting stress buildup is 
relieved through dislocation motion in the grains. There 
is a competition between the ordinary deformation mode 
(dislocations) and the grain boundary sliding. As the 
grain size is increased, the dislocation motion is eventu- 
ally expected to dominate, and we expect a transition to 
a behavior more like what is seen in coarse-grained ma- 
terials, incl. a normal Hall-Petch effect. The transition 
is beyond what can currently be simulated at the atomic 
scale, but we do see a weak increase in the dislocation 
activity when the grain size is increased: The increase in 
the fraction of hep atoms during a simulation is increas- 
ing shghtly with the grain size (Fig. O). 



VI. CONCLUSIONS 

Molecular dynamics and related techniques have been 
shown to be a useful approach to study the behavior of 
nanocrystalline metals. We have in detail investigated 
the plastic deformation of nanocrystalline copper, and 
shown that the main deformation mode is sliding in the 
grain boundaries. The sliding happens through a large 
number of small, apparently uncorrelated events, where 
a few grain boundary atoms (or a few tens of atoms) 
move past each other. It remains the main deformation 
mechanism at all grain sizes studied (up to 13 nm), even 
at zero temperature. As the grain boundaries are the 
main carriers of the deformation, decreasing the number 
of grain boundaries by increasing the grain size leads to 
a hardening of the material, a reverse Hall-Petch effect. 
This is observed in the simulations, both for T = OK and 
for T = 300K. 

The Young's moduli of the nanocrystalline systems are 
found to be reduced somewhat compared to the experi- 
mental value for polycrystalline copper with macroscopic 
grain sizes, decreasing with decreasing grain size. This 
indicates that the grain boundaries are elastically softer 
than the grain interiors. The Young's modulus is de- 
creasing with increasing temperature at a rate somewhat 
above what is seen experimentally in coarser-grained cop- 
per. 

Pores in the samples have a large effect on both the 
Young's modulus and the flow stress. This effect is en- 
hanced if the pores are mainly in the grain boundaries, 
as one could expect in samples produced experimentally 
by inert gas condensation. Sample porosity can explain a 
large number of experiments showing reverse Hall-Petch 
effect, but the softening due to grain boundary sliding 
may be important for high-quality samples with grain 
sizes close to the lower limit of what can be reached ex- 
perimentally. 
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FIG. 1. The initial configuration of a nanocrystalline cop- 
per system. The system contains approximately 100 000 
atoms arranged in 16 grains, giving an average grain diame- 
ter of 5.2 nm. Atoms are color-coded according to the local 
crystalline order. White atoms are in a perfect fee environ- 
ment. Green atoms are in local hep order, which for example 
corresponds to stacking faults. Atoms in any other environ- 
ment (in grain boundaries and dislocation cores) are colored 
purple. 
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FIG. 2. The radial distribution function g{r) for the system 
in Fig. |l| at T = 0. It is the average number of atoms per 
volume at the distance r from a given atom. The marks above 
the peaks indicate the positions of the (infinitely sharp) peaks 
in a perfect crystal at T = 0. The inset shows that g{r) does 
not vanish between the first and second peak as it does for 
a perfect crystal. This contribution comes from the grain 
boundaries. 
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FIG. 4. Stress-strain curves for nanocrystalline copper at 
K. Each curve is the average over seven simulations with 
a given average grain size, the error bars indicate the uncer- 
tainty of the average (Icr). Adapted from Ref. M. 
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FIG. 3. The grain size distributions obtained from the 
Voronoi construction. Grain size distributions are shown for 
simulations with intended average grain sizes of 3.28 and 
5.21 nm. In both cases the input configurations from four 
simulations were used to calculate the distributions. Fits 
to log-normal (solid) and normal (dashed) distributions are 
shown. Both distributions fit the data well. 
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FIG. 5. Stress-strain curves for nanocrystalline palladium 
at OK. Each curve is the average over two simulations with a 
given average grain size. 
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FIG. 6. Individual stress-strain curves for seven simula- 
tions of nanocrystalline copper with an average grain size of 
5.21 nm. 



FIG. 8. The final configuration of the nanocrystalline cop- 
per system shown in Fig. |l| after 10% deformation. Some 
stacking faults have appeared where partial dislocations have 
moved through the grains. The arrow marks a partial dislo- 
cation on its way through the grain 
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FIG. 7. Hall-Fetch plot of the simulations at K, showing 
a reverse Hall-Fetch relationship between the grain size and 
both the yield- and fiow stress. 
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FIG. 9. The radial distribution function for the system in 
Fig. y at T = 0. The marks above the peaks indicate the peak 
positions in a perfect fee crystal. 
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FIG. 10. The fraction of atoms in a given local crystal 
structure as a function of the grain size. Two curves are 
given for each type of atoms, one for the initial configura- 
tions (filled circles) and another showing how the fractions 
have changed after 10% deformation (open squares). We see 
how the fraction of atoms in the grain boundaries (marked 
"Other") decreases with increasing grain size, and how more 
atoms are in grain boundaries and stacking faults (hep) after 
the deformation than before. 
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FIG. 11. The same system as in Fig. H, but with the atoms 
colors according to their motion relative to the homogeneous 
deformation. Blue atoms have moved downwards relative to 
what would be expected if the deformation had been homoge- 
neous and elastic; red atoms have moved upwards. The scale 
indicate how far the atoms have moved (in Angstrom) dur- 
ing the 10% deformation. The grains are seen to be blue in 
the upper part, and red in the lower part, indicating that the 
grains have not been strained as much as the system. This in- 
dicates, that a major part of the deformation has been in the 
grain boundaries. Several dislocations have moved through 
the large grain in the upper-left corner. Their slip planes are 
clearly seen. 
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FIG. 12. Stress-strain curves for the system in Fig. hi, 
deformed at different temperatures. The strain rate (e) is 
5 X 10*s~^ except at OK, where no strain rate can be defined 
(see text). We see a clear softening at increasing tempera- 
tures. 
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FIG. 13. Stress-strain curves for nanocrystalline copper 
at 300 K for varying grain size. The curves for grain sizes 
d < 5.21 nm are the average of four simulations with a system 
size of 100 000 atoms and a strain rate of 5 x 10*s~^. The 
simulations with d > 6.56 nm were made with a system size 



of 1000 000 atoms and a strain rate of 10 s~ 

of this change in the strain rate is minimal, see section [V B 4 
and Fig. Pl9l One simulation was performed with d = 6.56 nm, 
two with d — 13.2 nm. The thermal fluctuations are less pro- 
nounced in the simulations with the larger system size. 
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FIG. 15. Yield and flow stresses as a function of tempera- 
ture. No yield stress is given for T = OK, as the yield stress 
at this temperature cannot dir ectly b e compared to that at 



The influence finite temperatures, see section 
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FIG. 14. Young's modulus as a function of temperature for 
a sample with an average grain size of 5.2 nm. 
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FIG. 16. A Hall-Fetch plot of the simulations at 300K. As 
at OK, a reverse Hall-Fetch relation is seen. 
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FIG. 17. The same system as in Fig. hll but after deforma- 
tion at 300 K, see the caption of that figure for a description 
of the color coding. The main deformation is seen to be in 
the grain boundaries, as was the case at OK. 
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FIG. 19. The effect of varying the strain rate. Stress-strain 
curves are shown for simulations of the same system at the 
same temperature (300 K), but with varying strain rates. A 
strong influence of strain rate is seen for strain rates above 



10« 



below that value the strain rate dependence is 



less pronounced. The simulations at the two lowest strain 
rates were stopped after 4% deformation. 
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FIG. 18. Development of the number of atoms in different 
local environments as a function of the strain, for a single 
simulation at OK and at 300 K. 
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FIG. 20. Summary of the effect of varying the strain rate. 
Both the yield stress and the flow stress are seen to vary with 
the strain rate (e), strongest for e > 10^s~^. 
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FIG. 21. Sample relaxation in the "elastic" region. 
The deformation is stopped after 0.4% deformation at 
e = 2.5 X 10''s~^ and T = 300ii'. The stress is seen to de- 
crease with a characteristic time of ~100 ps, stabilizing on 
a value of 2/3 the original stress. The thick black line is an 
exponential fit. 
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FIG. 22. Grain rotation as a function of strain and tem- 
perature. The rotation of five randomly selected grains were 
followed in simulations at different temperatures. The aver- 
age grain size is 5.2 nm. The error bars indicate the \a spread 
in the rotations. 
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FIG. 23. Effect of voids on the Young's Modulus (top) and 
on the flow stress (bottom) , showing a decrease of both quan- 
tities with increasing void fraction. The voids were generated 
by removing selected atoms, the "void fraction" is the frac- 
tion of atoms removed. Different methods were used to select 



the atoms to be removed, see section [VB(: 
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